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Abstract 

The present paper is devoted to the joint motion of two immiscible incompressible 
liquids in porous media. The liquids have different densities and initially separated by 
a surface of strong discontinuity (free boundary). We discuss the results of numerical 
simulations for exact free boundary problems on the microscopic level for the abso- 
lutely rigid solid skeleton and for the elastic solid skeleton of different geometries. The 
problems have a natural small parameter, which is the ratio of average pore size to the 
size of the domain in consideration. The formal limits as e \ results homogenized 
models, which are the Muskat problem in the case of the absolutely rigid solid skeleton, 
and the viscoelastic Muskat problem in the case of the elastic solid skeleton. The 
last model preserves a free boundary during the motion, while in the first model in- 
stead of the free boundary appears a mushy region, occupied by a mixture of two fluids. 
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§1. Introduction 

The RayleighTaylor instability, or RT instability (after Lord Rayleigh [1] and G. L Taylor 
[2]), is an instability of an interface between two fluids of different densities, which occurs 
when the lighter fluid is pushing the heavier fluid. The equivalent situation occurs when 
gravity is acting on two fluids of different density with the dense fluid above a fluid of lesser 
density. As the instability develops, downward-moving irregularities are quickly magnified 
into sets of inter-penetrating RayleighTaylor fingers. Therefore RT instability is sometimes 
qualified to be a fingering instability [3j. In the present paper we consider RT instability in 
an inhomogeneous fluid filling the pores in the solid skeleton. In addition to the undoubted 
theoretical interest, this problem is important for a number of practical problems. For 
example, the description of the displacement of one liquid by another in porous media. This 
problem still needs correct mathematical model. 

There are different types of mathematical models, but we are interested only in some of 
the fundamental models of continuum mechanics (such as, for example. Stokes equations for 
a slow motion of a viscous liquid, or Lame's equations for displacements of an elastic solid 
body), or in models asymptotically close to above mentioned ones. 
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Among mathematical models of a joint motion of two immiscible liquids the most 
trustable (or physically correct) one is the Muskat problem, suggested by M. Muskat[l]. 
This model describes a filtration of two immiscible incompressible liquids of different viscosi- 
ties and different densities, divided by some moving boundary (free boundary). The motion 
of the first liquid under the gravity in the domain f2+(t) with a constant viscosity and a 
constant density p'j is governed by the Darcy system of filtration 



k 



V 



VpJ + pJea, V-v+ = 0, xen+it), 



'1.1] 



for the macroscopic velocity and the pressure p~j of the liquid. The unit vector 62 coincides 
with the direction of the gravity. 

Correspondingly, the motion of the second liquid under the gravity in the domain Q^(t) 
with a constant viscosity fi and a constant density pj is governed by the Darcy system of 
filtration 

k 

V- = Vp]+pje2, V-v- = 0, xen-{t), (1.2) 

p 

for the macroscopic velocity v~ and the pressure . On the common free boundary r{t) = 
dQ^ (t) f] dQ^ (t) pressures and normal velocities are continuous: 



pj=Pf^ 



X e r(t), (1.3) 

■n = v- ■n = Vn, xeT{t), (1.4) 

where n is the unit normal vector to the boundary r{t) at the point 
G r(t) and Vn is the velocity in the normal direction of the boundary 
r(t) at the point x eT{t). 

Condition f ll.4p means that the boundary r(t) is a material surface. 
That is, it consists of the same set of material points during the motion. 
This fact permits a weak formulation of the Muskat problem. So, we 
define the pressure p/ of the inhomogeneous liquid as 

Pf = p^ if X e ^^^(t), Pf = pJ a X e ^^{t), 

the density p/ as 

Fig. 1. 1 - do- _ _ 

main Q+W, 2 - domain Pf = Pf if X G (t), Pf = pJ if X E (t), 

n- (t),3- free boundary g^^^ ^^le Vclocity V aS 

r(t). 

V = if X e i^^(t), V = V' if X e fi"(t). 
Then the unknown functions v, pf, and pf satisfy the Darcy system of filtration in the form 




V 



k 

--Vpf + Pf 62, V ■ V 
p 



0, X eVL, t > 0, 



and transport equations 



dpf _ dpf 
dt dt 



^ = ^ + 'U-Vp/ = 0, xeQ, t>0, 



;i.5) 



1.6) 
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The first equation in f ll.Sp (Darcy law) is understood in an usual sense almost everywhere in 
Qt = X (0, T), and the second equation (continuity equation) is understood in the sense of 
distribution. Transport equation is understood in a sense of distributions, if we use equalities 

V ■ V/i = V ■ (vp), V ■ V pf = V ■ (vpf). 

The problem is endowed with homogeneous boundary condition 

v-n = 0, xeS = dn, t>0, (1.7) 

where n is the normal vector to the boundary S, and initial conditions 

Pfix,0) = p'lix), xen, (1.8) 

with discontinuous initial data: 

p^f{x) = p^, if 03 G fi^ and p°(a;) = pj , if x e fi". 

So, one has two settings of the same Muskat problem. In both cases the problem is easy to 
formulate, but almost impossible to solve. For this reason, very little is known either about 
classical solutions or on weak solutions. There are only few results on a classical solvability 
locally in time or globally in time, but near explicit solutions, and there is no any result on 
a weak solvability (see, Ref. [5], Ref. [6], Ref. [7] and references there). 

Following R. Burridge and J. Keller Ref. [8], let us try to find more general physically 
correct mathematical models describing the same physical process. To explain idea we, at 
first, consider the Darcy system of filtration, which is responsible for the dynamics in the 
Muskat problem. It is well-known that this system is an asymptotic limit of the Stokes 
system for an incompressible viscous liquid, when dimensionless pore size tends to zero (see 
[8], [9j). R. Burridge and J. Keller have suggested to consider more general system 

= ^ ■ (x'«Mro(x, ^) + (1 - X')"A©(2;, W) - pi) + p' 62, (1.9) 

p + a.pV-w = {), (1.10) 

for the displacement w and pressure p of the continuum medium. The microscopic system 
01. 9p . fll.lOj) describes the joint motion of the viscous liquid in pore space and elastic solid 
skeleton and is understood in the sense of distributions. Roughly speaking, this system 
contains the Stokes system for the viscous liquid in the pore space the Lame's system 
for the solid skeleton in VLg and boundary condition (the continuity of the normal stresses) 
on the common boundary F = dVLj r\dVLa. In (11. 9p D(a;, w) is the symmetric part of Viu, 
is the characteristic function of the pore space £ = //L is the dimensionless pore size, 

_ L _ 2fx _ 2A 

gr^ rLgpo Lgpo 

p' = PfX' + P.(l - X'), ttp = c^IjX' + «J,.(1 - X'), 

I is an average size of pore, L is a characteristic size of the domain in consideration, r is a 
characteristic time of the process, pf and ps are the respective mean dimensionless densities 
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of the liquid in pores and the sohd skeleton correlated with the mean density of water po; 
apj and a^^s are the respective dimensionless speed of sound of the liquid in pores and the 
solid skeleton, g is the value of acceleration of gravity, /i is the viscosity of fluid, and A is the 
elastic Lame's constant. In what follows we suppose the structure of the domains VLj and 
f2s is periodic with a period e. That is x^i.'^) = with 1- periodic function xiv)- 

Theoretically the system (11.91) . f ll.lOp with corresponding boundary and initial conditions 
most accurately describes the given physical process, but it still has no any practical signif- 
icance. This one appears only after homogenization. So, we have to let all dimensionless 
criteria a,-, and ax to be variable functions, depending on the small parameter e, and 
find all limiting regimes as e — ?■ 0. 

It is clear that limiting regimes of the system (11. 9p . (ll.lOp depend on the dimensionless 
criteria 

ro = lima^, uo = \iman, Ao = limaA, 

e\0 e\0 s\0 

Cf = lim a„ f Cs = lim q;„ „, Ui = lim ^ 

(see [To]). For example, for filtration tq = 0, for absolutely rigid solid skeleton Xq = oo, and 
for incompressible media Cf = Cg = oo. Therefore, different sets of criteria lead to different 
homogenized models, asymptotically closed to the basic one. All of these models describe 
the same physical process, but with a different degree of approximation. 

For a filtration of an incompressible liquid ( c/ = oo) in an incompressible solid skeleton 
( Cs = oo) the standard homogenized model for fixed = /iq, oa = Aq, and ar = tq has a 
form 



^o^^72" = V ■ P + ^ 62, g = mpf + (1 - m)p^ 



dw /■* 
P = fio^i : ©(x, — ) + AoOla : ©(x, w) + m^it - r) : ©(x, w{x, T))dT. 
ot Jq 

But the same system ( 11. 9p . (ll.lOp on the microscopic level has another asymptotic limits, 
like Darcy system of filtration (tq = 0, C/ = Cs = oo, Aq = oo, po = 0, < pi < oo): 

_B/. _Vp + p^e2 , V - — =0, 
at Hi ot 

or Terzaghi - Biot system of poroelasticity (tq = 0, c/ = = oo, < Aq < oo, po = 0, 
< pi < oo): 

dw du 1 , „ 

^ ^dw . .du. 

V ■ (AoOIo : B{x, u)) -Vp + ge2 = 0, 
or the system of viscoelastic filtration (tq = 0, c/ = = oo, < Aq, Po < oo): 

V ■ P - Vp + ^ 62 = 0, V ■ w = 0, 
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dw r 

P = /ioOT4 : D(a;, — ) + AqOIs : D(x, w) + ^e{t - r) : D(x, 117(2;, r))cir 

Jo 

(see do], HU, [12]). 

The same scheme we may apply to the free boundary problem. On the microscopic level 
this model for filtration in incompressible media has the form 

V ■ {x'a^Bix, + (1 - x')axHx, w') - fl) + = 0, (1.11) 

V ■ = 0, (1.12) 



This problem has been completely studied in [13] . The main result there states, that for any 
e > there exists a weak solution {w'^, p^, p^} to the free boundary problem (11. lip - (I1.13P 
and for = /io, c^a = Aq the corresponding sequences converge as e \ to the solution 
{w, p, p} of the (homogenized) Muskat problem for the viscoelastic filtration 

V- P- Vj9 + pe9 = 0, V- = 0, ^ = 0. (1.14) 

dt 



The proof of this result has essentially used the notion of the two - scale convergence 

On the other hand the formal limit in (11.111) - (I1.13P as £ \ for Q!a — >■ 00 and = fiie'^ 
leads to the Muskat problem (II. 5p . (II. 6p . In is clear, that the same formal limit we obtain 
if as a basic model on the microscopic level we consider the free boundary problem for the 
Stokes system 

V ■ {pie'^3{x, - /I) + 62 = 0, xeQ},t> 0, (1.15) 

V ■ = 0, ^ = 0, xeQ},t>0, = 0,xe dQ}, t > 0, (1.16) 

only in the pore space fl'j. The last problem has been studied in |15j, where authors proved 
that for any e > the problem (11.150 has the unique classical solution {w'^, p^, p^}. 

The goal of the present paper is numerical simulations for the problem (11.110 - (I1.13P 
with = po, «A = Ao, and for the problem (I1.15p . (I1.16p . We have done numerical results 
for two different structures of the pore space: disconnected capillaries and disconnected solid 
skeleton in the unit square in M^. 

Numerical simulations of the problem in a single capillary in the absolutely rigid skeleton 
show the coincidence with results of [I6]. On the Figure 2 we may see the smooth free 
boundary (the surface of a strong discontinuity) in the capillary at different times. 
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t=0.003 f=0.255 t=0.525 t=1.035 
Fig. 2. The surface of a strong discontinuity at different times. 

The numerical results show that the motion of the liquids describing by fll.lSp . (I1.16P 
is affected by three factors: the ratio 6 = I of the densities and p~ of the liquids 
on the top and on the bottom, the viscosity p of fluids and the pore size e. Changing 
these parameters one gets different scenarios of Ray leigh- Taylor instability. For the problem 
(11 -lip - fll.lSp the process is affected by the same parameters 5, /i, e as before, and by the 
additional parameter A, which is the elastic Lame's constant. 

The limiting procedure {e \ 0) is modeled by increasing the number of capillaries for 
the first geometry, and the number of elementary squares for the second geometry. There- 
fore, we may assume that for sufficiently small e the problem f ll.lSp . fll.l6p describes the 
classical Muskat problem, while the problem (11. lip - (I1.13P describes the Muskat problem 
for viscoelastic filtration. 

Figure 3 shows numerical simulations for the first geometry for the model (ll.lSp . (I1.16P 
of the motion in an absolutely rigid solid skeleton (above), and for the model (II. lip - (I1.13P 
of the motion in the elastic solid skeleton (below). 




t=50 t=860 t=2631 t=3012 i=4873 

Fig. 3. Disconnected capillaries: numerical simulation for the absolutely rigid (above) and for the elastic solid skeleton. 
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t=5Q 



I ] 



t=860 



t=2631 



t=3012 



Fig. 4. The case of an elastic solid skeleton in a large scale. 

We also carried out numerical simulations for different values of A and 5: 

• for £ = 2 * 10~^, A \ and 5 — 1.25, there is a changing of the interface of fluids; 

• for £ = 2 * 10~^, A = 0.5 and 5\oo, there is a changing of the interface of fluids; 

• for 6 = 2* 10^^, A = 0.5 and 5 \ 1, there is no a changing of the interface of fluids. 

The same conclusion is valid for the second geometry (disconnected solid skeleton). Figure 
5 shows the comparative results for the same values S, /i, g, ps, A, e (see Table 1) and the 
same initial values. 




i=50 t=860 t=2631 t=3012 t=4873 

Fig. 5. Comparing results of absolutely rigid and elastic solid skeleton cases, absolutely rigid solid skeleton is above, elastic 

solid skeleton is down (disconnected skeleton). 



Parameter 


Absolutely rigid solid skeleton 


Elastic solid skeleton 




998.2 


998.2 


Pi 


800 


800 


Ps 




2000 


p+ 


10-2 


10-2 




9* 10-1 


9* 10-1 


A 




0.5 


e 


2 * 10-5 


2* 10-5 


L 


100 


100 



Table 1. The table of values for the cases of absolutely rigid solid skeleton and elastic solid skeleton. 



One may see, that for the same time of the process, in the Musket problem appears 
a mushy region (for definition see |17]), while the viscoelastic filtration preserves the free 
boundary. 

As for the first geometry, we have done numerical calculations for various values A and 

5: 

• for 5 = 1.25, e = 2 * 10^^ and A \ 0, there is a changing of the interface of fluids and 
the process of filtration is very slow; 

• for 6 \ oo, 6 = 2* 10-5 and A = 0.5, there is a changing of the interface of fluids and 
the process of filtration become faster with increasing 6; 

• for 6 = 1.25, e = 2* 10~^ and A \ oo, there is a changing of the interface of fluids and 
the process of filtration become faster with increasing A. 

On the basis of numerical calculations we can conclude that 
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1) for the liquid motion in the absolutely rigid solid skeleton instead of the free boundary 
appears a mushy region, while the motion in the elastic solid skeleton preserves the free 

boundary; 

2) when the liquids in the absolutely rigid solid skeleton completely change their positions, 
the liquids in the elastic solid skeleton still preserve their positions. 

Therefore, the solution to the model of a viscoelastic filtration is a classical one and 

possesses a smooth and stable free boundary, whereas the the solution to the Musket problem 
will at best be only a generalized solution with mushy region instead of the free boundary 



2 The problem statement 
2.1 Absolutely rigid skeleton 

We suppose that f2 is a unit square {0 < Xi < 1) x {0 < X2 < 1), Q = Q'^f U S"^ U Q^, where 
is a pore space (the domain occupied by the liquid), - the solid skeleton and S'^ is the 

"solid skeleton - pore space" interface. 

For the first geometry the pore space is a set of isolated capillaries 

n-l 

= IJ (^ek <xi < e{k + l)m) x (O < X2 <l). 

k=0 

Here e = 1/n, m is a porosity, < m < 1. 

For the second geometry the solid skeleton is a union of disjoint squares 

n— 1 

ni= [j{e{k + l)P<xi<e{k+l){l-P)) X {e{k + l)P <X2<e{k + l){l- P)), 
2/3 = 1- 

As we have mentioned above, two immiscible liquids are modeled by an inhomogeneous 
liquid, where the density p can take only two constant values p"*" or p~. The velocity v, 
the pressure pf and the density p/ of the inhomogeneous liquid in the pore space are 
described by the Stokes system 

V • (/xi£2D(a;, v) - pfl) + 62 = 0, V • v = 0, (2.1) 

and by the transport equation 

^ + v-Vpf = 0, (2.2) 

where 62 is the unit vector, which coincides with the direction of the gravity. 

Differential equations (2.1) and (2.2) are supplemented by the normalization condition 



/ pf{x,t)dx = 0, (2.3) 



the homogeneous boundary condition 

v = 0, xe SUS', t>0, (2.4) 
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where S = d^l, and the initial condition 

Pf{x,0)=p'}{x), xen}, (2.5) 
p'j(x) = p'f = const > 0, if 03 G = n 

and 

p^f{x) = pj = const > 0, a X E ^IJ = Vl'f n fi". 



2.2 Elastic solid skeleton 

The joint motion of the viscous hquid in pore space and elastic solid skeleton on the micro- 
scopic level is governed by the system fll.lOp - fll.lSp . which consists of the stationary Stokes 
system 

V ■ (/^oro(x,^)) - Vp/ + p/e2 = 0, V-^/ = 0, (2.6) 

for the displacements Wj and the pressure pj of the inhomogeneous liquid, the stationary 
Lame's equations 

V ■ (Ao B{w,)) - Vps + p. 62 = 0, V ■ 10, = 0, (2.7) 

for the displacements Wg and the pressure Ps of the elastic solid skeleton, and two boundary 
conditions on the common "solid skeleton - pore space" boundary S^: 

Wf = Ws, (2.8) 

duo 

{po ®(^) - P/I) ■ ^ = H-^s) - Psl) ■ n, (2.9) 

where n is the unit normal vector to the boundary S*^. 

The system (12. 6p - (12. 9p is supplemented with initial and boundary conditions 

Wf{x,0) = 0, xen'j, (2.10) 

w = 0, xeS, t>0, (2.11) 

where 

w{x^t) = Wf{x,t), ifx G ilf, t > 0, 

and 

w{x,t) = Ws{x,t), ifx G Qs, t > 0, 

and normalization condition 

/ p{x,t)dx = 0, (2.12) 

where 

p{x,t) = pf{x,t), ifa^Gfi/, t>0, 

and 

p{x, t) = Ps{x, t), if 33 G ^Is, t > 0, 

Next, we have to complete the Cauchy problem for the transport equation (12. 2p . On the 
microscopic level the transport equation for the liquid density is defined only in the pore 
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space For the exact microscopic model (see [10]) the characteristic function x"^ of the 
pore space is an unknown function and defined as a solution to the Cauchy problem 

It means that the solid - liquid interface 5*^ in the exact model is a material surface and 
we do not need boundary conditions for the liquid density on S"^. But our basic dynamic 
system here is a linear one, where = is the given function. Therefore the solid - liquid 
interface is no longer the material surface and we need the boundary condition for the 
liquid density on the part of S'^, where the liquid "enters" into the pore space. To avoid 
this, we extend the Cauchy problem ( 12. 2 p onto hole domain Q. At first, we suppose that the 
function p^j{x) is defined in hole domain f2 and 

p^f{x) = in fi^, p^f{x) = pj in Q~ . 

Finally, we rewrite the Cauchy problem for the liquid density as the Cauchy problem for the 
density of a mixture 

P = X'pf + (1 - X')ps 

in the form 

^ + ^-Vp = 0, {x,t)e Qr, p{x, 0) = p\x), xe Q, (2.13) 

where 

p^'\x) = x'p'f{x) + {l-xnPs. 



3 The computational algorithm 
3.1 The absolutely rigid solid skeleton 

As the numerical method for computer simulation fl2.ip - 02. 2p system of equations, sup- 
plemented with appropriate conditions fl2.3p - fl2.5p . we chose the large-particle method. 
This method is based on splitting the original differential equations in accordance with the 
physical processes they represent. The large-particle method is the development of Harlow's 
method of " particle- in-cell", which refers to the technique used to solve a certain class of 
partial differential equations, including the Rayleigh- Taylor instability for the Musket free 
boundary problem. 

The solution process for an evolution system (12. ip - (12. 2p is very difficult, because we 
have the presence of stratification. Thus the heterogeneity of the fluid requires additional 
calculation of the density field. To find the unknown functions the calculation process can 
be represented as a three-step scheme. 

Step 1: One has to calculate intermediate velocity v from the equation 

V ■ (pie^Biv)) + pe2 = 0, (3.1) 
and to find intermediate density p from the equation 
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^ + ^ • Vp = 0, (3.2) 

Step 2: The main difficulty in the numerical solution of equations (2.1), (2.2) associated 
with the calculation of the pressure field. The ffist significant success in overcoming these 
difficulties has been achieved through the idea of artificial compressibility [19j. It is very 
important to calculate unknown p from the equation: 

I + clfV ■v = 0. (3.3) 

Here we need to satisfy the solenoidal condition (V ■ v = 0). It is necessary to increase 
Cp each time step, while it will be equal to given accuracy. 

Step 3: We calculate the final values of the velocity v from the Stokes equation on the 
next time step: 

V ■ {fiie^B{v)) -Vp + pe2 = (3.4) 
and the density p from the equation: 

^ + v ■Vp = 0, (3.5) 

So, let us write the basic finite-difference scheme. The domain of integration is covered 
by a stationary (Eulerian) difference grid of arbitrary form (to abbreviate the exposition, a 
rectangular grid in a planar domain is considered, see Fig. 6): 



xj^+^Z^) =ihi, hi>0] i = 0, 1, A^i; 
= jh2, h2>0; j = 0,l,...,iV2; 



where hi, h2 is the size of the grid, A^^i, N2 are the numbers of grid cells, respectively, in 
the Xi and X2 direction (the point with coordinates matches with the center of the 

cell). Here, as in the original splitting method, we use the "checkerboard" grid. This makes 
it possible clearly interpret each cell as element of volume, which is characterized by a 
calculated pressure and density in it's center. 
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1' 



Xj.U 



Pl-U 

Pi-lj 



ptJ-1 



Pu 



Pli-i 
Pu-i 



,PI-I.i 
'Pi-lj 



J+I/2 



j-1/2 



i-3/2 hi i-1/2 



i+I/2 



j-3/2 



i+3/2 



Fig. 6. Template of the grid. 

Only quantities relating to the cell as a whole are varied, while the fluid is assumed to be 
momentarily constrained. For this reason, convective terms, corresponding to translational 
effects, are dropped from the equations. In the remaining equations, p is brought forward 
outside the differentiation symbol, and the equations (2.1)-(2.2) are solved for the derivatives 
of u, V with respect to the time. 

An elementary flnite-difference approximation of the equation (3.1) yields the following 
expressions: 



u 



i+3/2,j 



2 "r+l/2,j + ""-72 J 



"i+l/2,j+l 



9 7~/" _|_ 7~/" 

"'i+l/2,j ^ "1+1/2,^-1 



^i+l,j+l/2 



hi 

2 ^I!i+l/2 + ^l'-lj + l/2 



hi 



r,n _ O f.n _|_ ^.n 

-^i,j+3/2 ^ ^i,i+l/2 ^ ^i,j-l/2 Pi,j+l/2 



h{ /i2 /ii^^ 

where n - number of the time step, u, v - are intermediate values of the flow velocity. 
Quantities with fractional subscripts relate to cell boundaries, e.g.: 



'^i+l/2,j 



As in the large particles method, we need to calculate the flux density (3.2) through the 
cell boundaries. It is assumed throughout that the density of the large particle is in the 
motion only owing to the velocity component normal to the boundary. This values of the 
flow parameters in the next time layer are computed according to the following formulas 
(the direction of the flow is from left to right and upwards): 



21^1,] 



P: 



~\n+l 



i+l/2j 



T hi h2 

In cases, where we need to determine the function values in the grid points, not to meet 
their terms in Fig. 6, we used the arithmetic average, for example: 

_ 1^ ^ 
Pi+1/2 — ^yP-^+^d + Pi,j)- 
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At the final step (when a discrete model of the medium is being used) one should carry out 
an additional calculation of the density; this smoothness out fluctuations and increases the 
accuracy of the computations. Combining different representations of the steps, one obtains 
a series of difference schemes; this series, which constitutes the large-particle method. 

Chosen method may be interpreted from various points of view: the splitting method, 
the mixed EulerLagrange method, computation in local Lagrangian coordinates with scaling 
on the previous grid, the difference notation for conservation laws for a fluid element (large 
particle), and the Eulerian difference scheme. 

When we replace the differential problem by the flnite-difference representation, we should 
pay special attention to the approximation of boundary conditions, because their speciflc 
approximation affects the correctness of the method and stability of the scheme, as well as 
the velocity of convergence. 

The boundary conditions are formulated by introducing series of flctitious cells (so that 
every computation point becomes an interior point and the same algorithm is maintained 
for all cells). One layer is sufficient for the first-order approximation scheme, two layers are 
sufficient for the second order. As the result, in the case, where side walls are the solid 
surface, the impermeability condition is represented as 

ti_i/2,j = 0, (3.7) 

and the slip condition as 

Wi-i/2,j+i/2 = 0. (3.8) 

In the planar case, the geometrical characteristics of fractional cells may be determined 
by direct measurement. In the axially-symmetric case one needs an additional computation, 
incorporating the distance of each fractional cell from the axis of symmetry. The difference 
formulas for fractional cells are obtained by a slight modification of the difference formulas 
for whole cells. 

3.2 The elastic solid skeleton 

The computer simulation of joint motion of the viscous liquid (in the pore space) and elastic 
solid skeleton (at the microscopic level) is the numerical solution of the system, which consists 
of the stationary Stokes system fl2.6p for the displacement Wf and the pressure Pf of the 
inhomogeneous liquid and the stationary Lame's equations (12. 7p for the displacement Wg 
and the pressure ps of the elastic solid skeleton with appropriate conditions on the common 
border. 

As the numerical method to simulate the Rayleigh- Taylor instability with an elastic 
component of the skeleton, we chose the same method, which has been described above 
(large particles method). Considerate domain is replaced by the system of fluid particles, 
which match with the cell of Eulerian mesh (Fig. 6). 

The suggested algorithm consists of four main steps: 

Step 1. One has to solve the system of Lame's equations (12. 7p with given boundary and 
initial conditions for the displacement Wg and the solid pressure Ps'- 

Ws = Wf, cc G S"^, Wg = 0, X E S. (3.9) 
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The solution of this system is no different from the solution of the Stokes system (Sub- 
section 3.1). We enter the artificial compressibility and calculate unknown ps from the 
equation:: 

^ + cLV -^. = 0. (3.10) 

It is necessary to increase Cp^s each time step, while the solenoidal condition (V -Wg = 0) 
will be satisfied, since the fluid is incompressible. 

Also it should be noted, that the procedure for two-dimensional difference scheme of the 
Lame's equation is identical to the previous. And scheme itself is easier, unlike for Stokes 
equation, because we needn't to introduce the necessary differential approximation for the 
velocity. 

An elementary finite-difference approximation of the equation (2.7) yields the following 
expressions: 



w 



li+3/2,j 2 W^._^^/2j + '^li-/2,j _ '^li+l/2,j+l 2 ^i'l j+i/2J + W^ri+l/2j-l 



flf h 



w 



2i+lj+l/2 2 W2jj+i/2 + '"^2i-lJ+l/2 '"^2ij+3/2 ^ ""^2 j j+1/2 + ""^2 i,j-l/2 Psi,j+l/2 



where n - number of the time step, Wi, W2 are displacement values, ps is the density of the 
elastic solid skeleton. 

Step 2. Using known values Wg and Psi find normal stress on the boundary 5*^: 

(Ao©(ic«) - pj)n = A, (3.11) 
where n is the unit normal to the boundary S"^ . 



Step 3. Then solve the system of Stokes equations (12.61) in fi^ just as for absolutely rigid 
solid skeleton, repeating Step 1 - Step 3 from the previous subsection, with the condition on 
the common boundary S^: 

(/xoP(^) - Pfl)n = (XoH'^s) - Psl)n, (3.12) 

where we know the right side of the equality from the previous step. 

This boundary condition suggests, that the vector of displacement and pressure satisfies 
continuity of normal stresses on the common boundary between liquid and elastic solid 
skeleton. 

The found value of the fluid velocity replace to the transport equation ( 12.13^ in Q'j, 
where we find the density value pf for the next time step. 

Stage 4. When the velocity dwf/dt is known, in the liquid part we determine Wg on 
the next time step from the continuity of normal stresses and condition (2.8). 

Thus, considering behavior of fluids on the boundary of solid skeleton and, solving the 
system of equations ( 12. 6p . (12. 7p . ( I2.13P with appropriate initial and boundary conditions 
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(12. 8 p - fl2.12p . we obtain the numerical approximation of joint motion of fluid and elastic 
solid skeleton. 

Conclusions 

In the present paper we have shown, how to model physical processes using modern 
methods of the mathematical analysis. We started with the free boundary problem for a 
joint motion of two immiscible incompressible fluids on the microscopic level. Theoretically 
this mathematical model is the most suitable model, describing the given physical process. 
But this model has no practical value, because we have to solve the problem in the physical 
domain of some hundreds meters, while the coefficients oscillate on the physical size in some 
microns. The practical value of the model appears only after homogenization. In turn, the 
homogenization has at least two levels of approximation, which depend on the dimensionless 
criteria of the physical problem. The first level of approximation is the well - known Muskat 
problem. The second level of approximation of the free boundary problem on the microscopic 
level is the Muskat problem for a viscoelastic filtration. In our numerical calculations for the 
periodic structure, we simulated the homogenization by increasing the number of elementary 
cells per unit volume. The numerical results show that the solution to the Muskat problem 
is unstable, the free boundary louses its sharp structure and transforms into some domain 
(mushy region), occupied by the mixture of two liquids. Whereas, the solution to the Muskat 
problem for a viscoelastic filtration remains a classical one with a smooth free boundary. That 
is, the Muskat problem for a viscoelastic filtration is a natural generalization of the classical 
Muskat problem, which still remains unsolved as a mathematical task, and extremely difficult 
for the numerical realization. 

Acknowledgment 

This research is partially supported by the Federal Program "Research and scientific- 
pedagogical brainpower of Innovative Russia" for 2009-2013 (State Contract 02.740.11.0613). 

References 

[1] Rayleigh, Lord (John William Strutt). (1883) Investigation of the character of the equi- 
librium of an incompressible heavy fiuid of variable density. Proceedings of the London 
Mathematical Society 14: 170177. 

[2] Taylor, Sir Geoffrey Ingr.m (1950) The instability of liquid surfaces when accelerated in 
a direction perpendicular to their planes. Proceedings of the Royal Society of London. 
Series A, Mathematical and Physical Sciences 201 (1065): 192196. 

[3] H. B. Chen, B. Hilko, E.Panarella. (1994) The RayleighTaylorinstability in the spherical 
pinch Journal of Fusion Energy, V.13 (4): 275280. 

[4] M. Muskat. (1934) Two-fluid system in porous media. The encroachment of water into 
an oil sand. Physics 5 250 -264. 



16 



[5] Yi. Fahuai. (2003) Global classical solution of Muskat free boundary problem, J. Math. 
Anal. Appl. 288 442461. 



[6] E. Radkevich. (1995) On the spectrum of the pencil in the Verigin- Muskat problem. 
Sbornik: Mathematics 80 No.l 33- 74. 

[7] M.Siegel, R.E.Caflish, S.Howison. (2004) Global existence, singular solutions, and ill- 
posedness for the Muskat problem. Comm. on Pure and Appl. Math. LVII 1 - 38. 

[8] R. Burridge and J. B. Keller. (1981) Poro elasticity equations derived from microstructure. 
Journal of Acoustic Society of America 70 No. 4 1140 - 1146. 

[9] E. Sanchez-Palencia. (1980) Non-Homogeneous Media and Vibration Theory, Lecture 
Notes in Physics, Vol. 129, (Springer- Verlag). 

[10] A. Meirmanov. (2007) Nguetseng's two-scale convergence method for filtration and seis- 
mic acoustic problems in elastic porous media, Siberian Mathematical Journal 48, 519 - 
538. 

[11] K. Terzaghi. (1923) Die Berechnung der Durchlassigkeitsziffer des Tones aus dem Ver- 
lauf der hydrodynamischen Spannungsercheinungen, Sitzung berichte. Akademie der Wis- 
senschaften, Wien Mathematiesch-Naturwissenschaftliche Klasse 132 No.IIa 104 - 124. 

[12] M. Biot. (1941) General theory of three dimensional consolidation. Journal of Applied 
Physics 12 155 - 164. 

[13] A. Meirmanov. The Muskat problem for a viscoelastic filtration. Accepted for publication 
to Interfaces and Free Boundaries. 

[14] G. Nguetseng. (1989) A general convergence result for a functional related to the theory 
of homogenization. SIAM J. Math. Anal. 20, 608 - 623. 

[15] S. Antontsev, A. Meirmanov, B. Yurinsky. (2000) A Free Boundary Problem for Stokes 
Equations: Glassical Solution. Interfaces and Free Boundaries 2, 413-424. 

[16] B. J. Daly. (1967) Numerical study of two fiuid Raylaigh- Taylor instability. - Phys. 
Fluids, No.2, 297-307. 

[17] A. Meirmanov. (1992) The Stefan problem, Walter de Gruyter, Berlin-New York, 244pp. 

[18] O. M. Bclotscrkovski. (1984) Numerical simulation in continuum mechanics.- Moscow, 
Science, 510-520. 

[19] N. N. Vladimirova, B. G. Kuznetsov, N. N. Yanenko. (1996) Numerical calculations 

of the symmetric fiow around a fiat plate by the viscous incompressible fiuid. - Some 
problems in Computational and Applied Mathematics. Novosibirsk, 186-192. 



17 



